EBAIIn2 - ChIP-seq Workshop
Elodie Darbo\(^{1}\) - elodie.darbo@u-bordeaux.fr; Stephanie Le Gras\(^{2}\) - slegras@igbmc.fr
2021 May
\(^1\) Université de Bordeaux, INSERM U1218
\(^2\) GenomEast platform, IGBMC
1 Introduction
During this training session, we are going to analyze data from Laurette et al. (Laurette et al. 2015) which data are deposited in GEO as GSE61967. There are ChIP-seq and RNA-seq data. Data were downloaded from GEO and processed.
1.1 Define the paths to the data
We create a list containing the absolute paths to the data we will use along the practice and their file names.
params <- list(
chipseq= "/shared/projects/ebai2021_n2/data/chip_seq/ChIPseq/",
mapping.chipseq= "stats_mappingChIPseq.tsv",
peakcalling.chipseq="stats_peakCalling.tsv",
brg1= "BRG1siCTRL_CHIP-seq_peaks.narrowPeak",
mitf= "MITF_CHIP-seq_peaks.narrowPeak",
sox10= "SOX10_CHIP-seq_peaks.narrowPeak",
brg1.bw= "BRG1siCTRL_CHIP-seq.bw",
mitf.bw= "MITF_CHIP-seq.bw",
sox10.bw= "SOX10_CHIP-seq.bw",
rnaseq= "/shared/projects/ebai2021_n2/data/chip_seq/RNAseq/RNAseq_diff_norm.RData"
)1.2 Processing of RNA-seq data
1.2.1 Preprocessing
Reads were preprocessed in order to remove adapter, polyA and low-quality sequences (Phred quality score below 20). After this preprocessing, reads shorter than 40 bases were discarded for further analysis. These preprocessing steps were performed using cutadapt (Martin 2011) version 1.10.
1.2.2 Mapping
Reads were mapped onto the hg38 assembly of Homo sapiens genome using STAR (Dobin et al. 2013) version 2.5.3a.
1.2.3 Quantification
Gene expression quantification was performed from uniquely aligned reads using htseq-count verson 0.6.1.p1 (Anders, Pyl, and Huber 2015), with annotations from Ensembl version 103 and “union” mode. Only non-ambiguously assigned reads have been retained for further analyses.
1.2.4 Normalization
Read counts have been normalized across samples with the median-of-ratios method proposed by Anders and Huber (Anders and Huber 2010), to make these counts comparable between samples.
1.2.5 Differential expression analysis
Differential expressed genes between wild type ans shBRG1 cells was performed using the edgeR R package (Robinson, McCarthy, and Smyth 2010). We called significant changes when FDR < 0.01, absolute log fold change over 1 and minimum average log normalized count over 5.
1.3 Processing of ChIP-seq data
1.3.1 Mapping
Reads were mapped to Homo sapiens genome (assembly hg38) using Bowtie (Langmead et al. 2009) v1.0.0 with default parameters except for “-p 3 -m 1 –strata –best –chunkmbs 128”. The following table shows the number of reads aligned to the Homo sapiens genome.
| Sample ID | Sample name | raw reads | aligned reads | multimapped reads | unmapped reads |
|---|---|---|---|---|---|
| SRR1594290 | MITF_CHIP-seq | 84856252 | 54369181 (64.07%) | 12823841 (15.11%) | 17663230 (20.82%) |
| SRR1594291 | SOX10_CHIP-seq | 38202152 | 29594694 (77.47%) | 6128890 (16.04%) | 2478568 (6.49%) |
| SRR1594292 | BRG1siCTRL_CHIP-seq | 50190992 | 40134919 (79.96%) | 6403809 (12.76%) | 3652264 (7.28%) |
| SRR1594293 | BRG1siMITF_CHIP-seq | 64765675 | 44480241 (68.68%) | 9357361 (14.45%) | 10928073 (16.87%) |
| SRR1594294 | GFPsiCTRL_CHIP-seq | 38148419 | 26273410 (68.87%) | 5722823 (15.00%) | 6152186 (16.13%) |
| SRR1594295 | MITF_Input | 29433042 | 19970925 (67.85%) | 4049711 (13.76%) | 5412406 (18.39%) |
| SRR1594296 | SOX10_Input | 35449561 | 27381422 (77.24%) | 7045192 (19.87%) | 1022947 (2.89%) |
| SRR1596499 | BRG1siSOX10_CHIP-seq | 42745544 | 34418403 (80.52%) | 6855741 (16.04%) | 1471400 (3.44%) |
1.3.2 Peak Calling
Prior to peak calling, reads falling into Encode blacklisted regions (“(2014) Mod/Mouse/humanENCODE: Blacklisted Genomic Regions for Functional Genomics Analysis - Anshul Kundaje” n.d.) were removed using bedtools intersect v2.26.0 (Quinlan and Hall 2010). Then peak calling was done with Macs2 v2.1.1 with default parameters.
| IP sample | Input sample | No. of peaks |
|---|---|---|
| BRG1siCTRL_CHIP-seq | GFPsiCTRL_CHIP-seq | 72024 |
| BRG1siMITF_CHIP-seq | GFPsiCTRL_CHIP-seq | 33984 |
| BRG1siSOX10_CHIP-seq | GFPsiCTRL_CHIP-seq | 73267 |
| MITF_CHIP-seq | MITF_Input | 9702 |
| SOX10_CHIP-seq | SOX10_Input | 4538 |
1.3.3 Generation of BigWig files
Normalized BigWig files were generated using Homer (Heinz et al. 2010) makeUCSCfile v4.11.0 with the following parameter ’-norm 1e7’ meaning that data were normalized to 10M reads.
2 Description and statistics of BRG1 dataset
2.1 Load data
Peak files are in narrowPeak format which is of the form (source):
- chrom - Name of the chromosome (or contig, scaffold, etc.).
- chromStart - The starting position of the feature in the chromosome or scaffold. The first base in a chromosome is numbered 0.
- chromEnd - The ending position of the feature in the chromosome or scaffold. The chromEnd base is not included in the display of the feature. For example, the first 100 bases of a chromosome are defined as chromStart=0, chromEnd=100, and span the bases numbered 0-99.
- name - Name given to a region (preferably unique). Use “.” if no name is assigned.
- score - Indicates how dark the peak will be displayed in the browser (0-1000). If all scores were "‘0"’ when the data were submitted to the DCC, the DCC assigned scores 1-1000 based on signal value. Ideally the average signalValue per base spread is between 100-1000.
- strand - +/- to denote strand or orientation (whenever applicable). Use “.” if no orientation is assigned.
- signalValue - Measurement of overall (usually, average) enrichment for the region.
- pValue - Measurement of statistical significance (-log10). Use -1 if no pValue is assigned.
- qValue - Measurement of statistical significance using false discovery rate (-log10). Use -1 if no qValue is assigned.
- peak - Point-source called for this peak; 0-based offset from chromStart. Use -1 if no point-source called.
## The package ChIPseeker provides a function to load peak files such as narrowPeaks as GRanges objects
## Here BRG1 peak set is loaded into a list of peaks
## this list can be extended if there are more datasets
library(ChIPseeker)
peaks <- list()
peaks[["BRG1"]] <- readPeakFile(paste0(params$chipseq, params$brg1), as="GRanges")
peaks## $BRG1
## GRanges object with 72024 ranges and 7 metadata columns:
## seqnames ranges strand | V4 V5
## <Rle> <IRanges> <Rle> | <character> <integer>
## [1] chr1 980402-981178 * | BRG1siCTRL_CHIP-seq_.. 77
## [2] chr1 983476-984526 * | BRG1siCTRL_CHIP-seq_.. 108
## [3] chr1 1000729-1001179 * | BRG1siCTRL_CHIP-seq_.. 90
## [4] chr1 1001762-1002054 * | BRG1siCTRL_CHIP-seq_.. 40
## [5] chr1 1020914-1021207 * | BRG1siCTRL_CHIP-seq_.. 114
## ... ... ... ... . ... ...
## [72020] chrY 18936757-18937116 * | BRG1siCTRL_CHIP-seq_.. 37
## [72021] chrY 19577774-19578316 * | BRG1siCTRL_CHIP-seq_.. 73
## [72022] chrY 19891214-19892094 * | BRG1siCTRL_CHIP-seq_.. 114
## [72023] chrY 19892553-19893192 * | BRG1siCTRL_CHIP-seq_.. 55
## [72024] chrY 21837647-21838039 * | BRG1siCTRL_CHIP-seq_.. 55
## V6 V7 V8 V9 V10
## <character> <numeric> <numeric> <numeric> <integer>
## [1] . 5.81321 10.22971 7.71662 166
## [2] . 5.97792 13.59661 10.83201 384
## [3] . 6.65677 11.62894 9.01201 199
## [4] . 3.78100 6.17367 4.02109 209
## [5] . 8.11634 14.29797 11.46427 158
## ... ... ... ... ... ...
## [72020] . 4.53175 5.93549 3.79936 167
## [72021] . 6.32404 9.87728 7.38650 165
## [72022] . 7.86030 14.29797 11.46427 711
## [72023] . 5.29987 7.83781 5.50942 390
## [72024] . 5.55592 7.83781 5.50942 152
## -------
## seqinfo: 24 sequences from an unspecified genome; no seqlengthsPeaks are stored as GenomicRanges objects; this is an R format which look like the bed format, but is optimized in terms of memory requirements and speed of execution.
We can start by computing some basic statistics on the peak sets.
2.2 How many peaks were called?
Compute the number of peaks per dataset.
# sapply() function takes list, vector or data frame as input and gives output in vector or matrix
# sapply apply the same function (here length) to all elements
# of the list "peaks"
sapply(peaks,length)## BRG1
## 72024Make a simple barplot showing the number of peaks per chipped protein.
barplot(sapply(peaks,length))
Let’s create a barplot with ggplot2 out of this data.
# Load ggplot2 library
library(ggplot2)
# create a table with the data to display
peak.lengths <- data.frame(IP=names(peaks),
NbPeaks=sapply(peaks,length))
# make the barplot
ggplot(peak.lengths, aes(x=IP, y=NbPeaks)) +
geom_bar(stat="identity")
# Let's add colors to the barplot
# In R it exists some already defined colors palettes
# the most widely used palette is RColorBrewer.
# This R library offers several color palettes
# See:
library(RColorBrewer)
par(mar=c(3,4,2,2))
display.brewer.all()
# now lets add colors to the barplot
# first ass the new information, fill=IP to let ggplot know
# that colors change based on chipped protein
ggplot(peak.lengths, aes(x=IP, y=NbPeaks, fill=IP)) +
geom_bar(stat="identity")
# if we want to use colors from RColorBrewer library
# with the "Set1" color palette
ggplot(peak.lengths, aes(x=IP, y=NbPeaks, fill=IP)) +
geom_bar(stat="identity")+
scale_fill_brewer(palette="Set1")
2.3 How large are these peaks?
Compute statistics on BRG1 peak sizes
## we use the function width() from GenomicRanges
library(GenomicRanges)
summary(width(peaks$BRG1))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 254.0 308.0 427.0 578.7 671.0 9186.0Create a simple boxplot of the peak sizes.
peak.width <- lapply(peaks,width)
boxplot(peak.width)
Now, create a nice looking boxplot with ggplot2. ggplot takes a data frame as input. We can either create a date frame, this is what we have done when we created the barplot. Here, we are going to use a package that creates a data frame from other types of data: reshape2.
# Load the package
library(reshape2)
peak.width.table <- melt(peak.width)
head(peak.width.table)## value L1
## 1 777 BRG1
## 2 1051 BRG1
## 3 451 BRG1
## 4 293 BRG1
## 5 294 BRG1
## 6 351 BRG1## create boxplot
ggplot(peak.width.table, aes(x=L1, y=value)) +
geom_boxplot()
Enhance it.
# - theme_classic() change grey background to white background
# - fill=L1 and scale_fill_brewer(palette="Set1") colors boxplots
# based on chipped protein and with colors from RColorBrewer Set1 palette
# - labs changes x and y axis labels and legend title
# - scale_y_log10() set y axis to a log scale so that we can have a nice
# view of the data in small values
ggplot(peak.width.table, aes(x=L1, y=value, fill=L1)) +
geom_boxplot()+
theme_classic()+
scale_fill_brewer(palette="Set1")+
labs(x = "", y = "Peak sizes", fill = "")+
scale_y_log10()
2.4 Peak filtering
To make sure we keep only high quality data. We are going to select one those peaks having a qValue >= 8. The qValue corresponds to the 9th column of narrowPeak files. So, we are going to set a threshold on this.
## Select high quality peaks
peaks$BRG1 <- peaks$BRG1[peaks$BRG1$V9 >= 8,]
## Compute the number of remaining peaks
length(peaks$BRG1)## [1] 308742.5 Where are the peaks located over the whole genome?
Sometime, peaks may occur more in some chromosoms than others. We can display the genomic distribution of the peaks along the chromosomes, using the covplot function from ChIPSeeker. Height of peaks is drawn based on the peak scores.
# genome wide BRG1 peak distribution
covplot(peaks$BRG1, weightCol="V5")
# chromosome wide BRG1 peak distribution
covplot(peaks$BRG1, chrs=c("chr1", "chr2"), weightCol="V5")
2.6 Functional annotation: genomic features enriched in BRG1 peaks
2.6.1 Generate annotation
We can assign peaks to the closest genes and genomic features (introns, exons, promoters, distal regions, etc…). This is done using the function annotatePeak which compares peak positions with the genomic feature positions of the reference genome. This function returns a complex object which contains all this information.
## org.Hs.eg.db is an R object that contains mappings between Entrez Gene identifiers and GenBank accession numbers.
library(org.Hs.eg.db)
## Annotate peaks for all datasets and store it in a list
## Here TSS regions are regions -1000Kb/+100b arount TSS positions
## Peak annotations are stored in a list
## Load transcript annotation
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
peakAnno <- list()
peakAnno[["BRG1"]] = annotatePeak(peaks$BRG1, tssRegion=c(-1000, 100), TxDb=txdb, annoDb="org.Hs.eg.db")## >> preparing features information... 2021-05-21 19:50:19
## >> identifying nearest features... 2021-05-21 19:50:21
## >> calculating distance from peak to TSS... 2021-05-21 19:50:22
## >> assigning genomic annotation... 2021-05-21 19:50:22
## >> adding gene annotation... 2021-05-21 19:50:57
## >> assigning chromosome lengths 2021-05-21 19:50:57
## >> done... 2021-05-21 19:50:57## Visualize and export annotation as a data table
# as.data.frame(peakAnno$BRG1)
head(as.data.frame(peakAnno$BRG1))## seqnames start end width strand V4 V5 V6
## 1 chr1 983476 984526 1051 * BRG1siCTRL_CHIP-seq_peak_2 108 .
## 2 chr1 1000729 1001179 451 * BRG1siCTRL_CHIP-seq_peak_3 90 .
## 3 chr1 1020914 1021207 294 * BRG1siCTRL_CHIP-seq_peak_5 114 .
## 4 chr1 1069156 1070192 1037 * BRG1siCTRL_CHIP-seq_peak_10 96 .
## 5 chr1 1079517 1080087 571 * BRG1siCTRL_CHIP-seq_peak_11 179 .
## 6 chr1 1304853 1306013 1161 * BRG1siCTRL_CHIP-seq_peak_15 108 .
## V7 V8 V9 V10
## 1 5.97792 13.59661 10.83201 384
## 2 6.65677 11.62894 9.01201 199
## 3 8.11634 14.29797 11.46427 158
## 4 6.90466 12.34852 9.65604 410
## 5 9.95583 21.25171 17.96153 213
## 6 5.97792 13.59661 10.83201 176
## annotation geneChr geneStart geneEnd
## 1 Distal Intergenic 1 975204 982093
## 2 Promoter 1 1001138 1014540
## 3 Intron (ENST00000379370.7/375790, intron 1 of 35) 1 1020123 1056118
## 4 Exon (ENST00000412397.2/100288175, exon 10 of 10) 1 1070967 1074306
## 5 Downstream (1-2kb) 1 1082146 1084072
## 6 Promoter 1 1292390 1305929
## geneLength geneStrand geneId transcriptId distanceToTSS ENSEMBL
## 1 6890 2 84808 ENST00000341290.6 -1383 ENSG00000187642
## 2 13403 1 9636 ENST00000624697.4 0 ENSG00000187608
## 3 35996 1 375790 ENST00000620552.4 791 ENSG00000188157
## 4 3340 2 401934 ENST00000453464.3 4114 ENSG00000237330
## 5 1927 2 54991 ENST00000464905.1 3985 ENSG00000131591
## 6 13540 2 116983 ENST00000492936.5 0 ENSG00000131584
## SYMBOL GENENAME
## 1 PERM1 PPARGC1 and ESRR induced regulator, muscle 1
## 2 ISG15 ISG15 ubiquitin like modifier
## 3 AGRN agrin
## 4 RNF223 ring finger protein 223
## 5 C1orf159 chromosome 1 open reading frame 159
## 6 ACAP3 ArfGAP with coiled-coil, ankyrin repeat and PH domains 3All the peak information contained in the peak list will be retained in the output of annotatePeak. The position and strand information of nearest genes are reported. The distance from peak to the TSS of its nearest gene is also reported. The genomic region of the peak is reported in annotation column. Since some annotation may overlap, ChIPseeker adopted the following priority in genomic annotation :
- Promoter
- 5’ UTR
- 3’ UTR
- Exon
- Intron
- Downstream
- Intergenic
- Downstream is defined as the downstream of gene end.
This hierachy can be customized using the parameter genomicAnnotationPriority.
annotatePeak report detail information when the annotation is Exon or Intron, for instance “Exon (uc002sbe.3/9736, exon 69 of 80)”, means that the peak is overlap with an Exon of transcript uc002sbe.3, and the corresponding Entrez gene ID is 9736 (Transcripts that belong to the same gene ID may differ in splice events), and this overlaped exon is the 69th exon of the 80 exons that this transcript uc002sbe.3 prossess.
Parameter annoDb is optional, if provided, extra columns including SYMBOL, GENENAME, ENSEMBL/ENTREZID will be added. The geneId column in annotation output will be consistent with the geneID in TxDb. If it is ENTREZID, ENSEMBL will be added if annoDb is provided, while if it is ENSEMBL ID, ENTREZID will be added.
–
Reminder: The TxDb class is a container for storing transcript annotations.
Bioconductor provides several packages containing TxDb objects for model organisms sur as Human and mouse. For instance, TxDb.Hsapiens.UCSC.hg38.knownGene, TxDb.Hsapiens.UCSC.hg19.knownGene for human genome hg38 and hg19, TxDb.Mmusculus.UCSC.mm10.knownGene and TxDb.Mmusculus.UCSC.mm9.knownGene for mouse genome mm10 and mm9, etc.
User can also prepare their own TxDb by retrieving information from UCSC Genome Bioinformatics and BioMart data resources by R function makeTxDbFromBiomart and makeTxDbFromUCSC.
One can also create a TxDb objects for his favourite organism using an annotation file in GTF/GFF format using the function makeTxDbFromGFF or the package GenomicFeatures.
Expand to find Coturnix japonica example
## download GTF file
download.file("https://ftp.ncbi.nlm.nih.gov/genomes/all/annotation_releases/93934/101/GCF_001577835.2_Coturnix_japonica_2.1/GCF_001577835.2_Coturnix_japonica_2.1_genomic.gtf.gz", "Coturnix_japonica_2.1.annotation.gtf.gz")
## Build TxDb object
library(GenomicFeatures)
txdb = makeTxDbFromGFF("Coturnix_japonica_2.1.annotation.gtf.gz")
## To save the txdb database
library(AnnotationDbi)
saveDb(txdb, 'txdb.Coturnix_japonica_2.1.sqlite')
## load it when needed
library(AnnotationDbi)
txdb = loadDb(file = 'txdb.Coturnix_japonica_2.1.sqlite')–
2.6.2 Visualize genomic annotation
We can now analyze more in details genomic features associated to our peaks (introns, exons, promoters, distal regions,…).
## distribution of genomic features for BRG1 peaks
# as a pie chart - which is the most widely used representation in publication
plotAnnoPie(peakAnno$BRG1)
# as a barplot
plotAnnoBar(peakAnno$BRG1)
However since some annotation overlap, ChIPseeker provides functions that help having a view of full annotation overlap.
library(UpSetR)
upsetplot(peakAnno$BRG1)
2.6.3 Importance of the annotation database
In the previous section, we have annotated peaks using annotation from UCSC’s knownGene. The database chosen for annotation can have an impact on subsequent results including data integration such as comparison between genes associated to peaks and gene expression using RNA-seq data for instance.
Let’s look at the overlap between genes expressed in RNA-seq data with genes associated to peaks (let’s remind that peaks are associated to closest genes if no other evidences are used).
## Download GTF files
library(GenomicFeatures)
## GTF file downloaded from https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_full_analysis_set.refseq_annotation.gtf.gz
txdb.refseq = makeTxDbFromGFF("data/GCA_000001405.15_GRCh38_full_analysis_set.refseq_annotation.gtf.gz")
txdb.ensembl = makeTxDbFromGFF("data/Homo_sapiens.GRCh38.103_UCSC_chr.gtf.gz")
## Ensembl TxDB can also been download using the function makeTxDbFromBiomart
# txdb.ensembl <- makeTxDbFromBiomart(biomart="ENSEMBL_MART_ENSEMBL",
# host="nov2020.archive.ensembl.org", dataset="hsapiens_gene_ensembl")
annot <- list()
annot[["ensembl"]] <- annotatePeak(peaks$BRG1, tssRegion=c(-1000, 100), TxDb=txdb.ensembl, annoDb="org.Hs.eg.db")
annot[["refseq"]] <- annotatePeak(peaks$BRG1, tssRegion=c(-1000, 100), TxDb=txdb.refseq, annoDb="org.Hs.eg.db")
rnaseq <- read.table(paste0(params$path.rnaseq, params$rnaseq), header=TRUE,
check.names = FALSE, quote="", sep="\t")2.7 Heatmaps: visualisation of binding profiles
Heatmap are widely used representation of ChIP-seq data as it allows visualization of read enrichment at various locations at the same time. For instance, one may want to represent read of its chipped protein enrichment +/-5Kb around all TSS of the reference genome or compare read enrichment at the same locations in many chip-seq datasets.
2.7.1 BRG1 peak distribution
We want to know if BRG1 is binding large and/or narrow regions, unique and/or tandem etc.
First of all, for calculating the profile of ChIP peaks binding to the BRG1 center, we need to define peak centers and extend them each side. Then align the peaks that are mapped to these regions, and generate the tagMatrix.
## Load the library
library(GenomicRanges)
# Compute the center of the peaks and attribute it to a new
# column in the metadata of the BRG1 peak GRanges
peaks$BRG1$center.peak <- (start(peaks$BRG1) + end(peaks$BRG1))/2
# For computation and memory efficiency reasons,
# we subset the top 10K peaks according to the FDR column (V9)
top.10000 <- peaks$BRG1[order(peaks$BRG1$V9,decreasing=T)][1:10000]
# Generate peak center GRanges for the 10K top peaks
centers.BRG1 <- GRanges(seqnames(top.10000),
IRanges(start=top.10000$center.peak,
end=top.10000$center.peak))
# Extend each side of 2000 bp
extended.2K.BRG1 <- centers.BRG1
ranges(extended.2K.BRG1) <- IRanges(start=start(centers.BRG1)-2000,
end=end(centers.BRG1)+2000)
## compute the density of peaks within the promoter regions
tagMatrix <- getTagMatrix(peaks$BRG1, windows=extended.2K.BRG1)
## plot the density
tagHeatmap(tagMatrix, xlim=c(-2000, 2000), color="red")
The regions are ordered relatively to their peak enrichment.
We can visualize a summary of the binding profiles looking at the corresponding average profile. This kind of profiles is much less greedy, we can extend a bit more (e.g. +/- 5000) from the peak centers redoing the previous steps. Definition of the regions have to be redone. Try to create the GRanges and the tagMatrix.
“Click to expand and see the code”
# Extend from the peak center
extended.5K.BRG1 <- centers.BRG1
ranges(extended.5K.BRG1) <- IRanges(start=start(centers.BRG1)-5000,
end=end(centers.BRG1)+5000)
## compute the density of peaks within the promoter regions
tagMatrix <- getTagMatrix(peaks$BRG1, windows=extended.5K.BRG1)
# Plot the profile
plotAvgProf(tagMatrix, xlim=c(-5000, 5000),
xlab="Distance to peak center", ylab = "Peak Count Frequency")
It looks like BRG1 is having several binding patterns but the binary nature of the signal (presence/absence of peaks) and the ordering of the rows do not allow us to appreciate them.
2.7.2 Read enrichment in BRG1 peaks
For computation and memory efficiency reasons, we are not going to look at read coverage of each position in the regions of interest but they are rather split into non-overlapping equally sized windows. Thus, we need to build a matrix composed of rows that are all BRG1 peaks and columns that contain read enrichment in all bins.
This allows to consider a bigger set of peaks and covered region. We will now analyze the whole set of BRG1 peaks over 10Kb (+/- 5000bp). We will ask for a hundred bins each side of the center resulting in 200 windows of 50 bp.
2.7.2.1 Prepare the signal matrices
We need to load bigwig files for all datasets that we want to visualize. Data are imported using a function from the rtracklayer package.
# # Generate peak center GRanges for all the peaks
centers.BRG1 <- GRanges(seqnames(peaks$BRG1),
IRanges(start=peaks$BRG1$center.peak,
end=peaks$BRG1$center.peak))
# Extend each side of 5000 bp
extended.5K.BRG1 <- centers.BRG1
ranges(extended.5K.BRG1) <- IRanges(start=start(centers.BRG1)-5000,
end=end(centers.BRG1)+5000)
# load the library
library(rtracklayer)
# load the bw file for BRG1
cvg.BRG1 <- import(file.path(params$chipseq,params$brg1.bw),format="BigWig",
which=extended.5K.BRG1,
as="RleList")
cvg.BRG1## RleList of length 25
## $chr1
## numeric-Rle of length 248956422 with 1333289 runs
## Lengths: 979000 3 48 15 32 ... 1288 73 36 22594
## Values : 0.00 0.32 0.16 0.49 0.32 ... 0.00 0.17 0.33 0.00
##
## $chr10
## numeric-Rle of length 133797422 with 510080 runs
## Lengths: 368603 154 85 154 304 ... 90 58 154 272220
## Values : 0.00 0.16 0.00 0.16 0.00 ... 0.65 0.00 0.33 0.00
##
## $chr11
## numeric-Rle of length 135086622 with 542307 runs
## Lengths: 202802 21 117 221 41 ... 9 70 154 74826
## Values : 0.00 0.49 0.16 0.00 0.16 ... 0.16 0.00 0.16 0.00
##
## $chr12
## numeric-Rle of length 133275309 with 564715 runs
## Lengths: 194505 18 115 154 136 ... 22 68 19 189737
## Values : 0.00 0.32 0.00 0.16 0.00 ... 0.16 0.65 0.49 0.00
##
## $chr13
## numeric-Rle of length 114364328 with 338291 runs
## Lengths: 19776001 19 5 12 ... 10 13 43427
## Values : 0.00 1.29 1.13 0.97 ... 0.16 0.65 0.00
##
## ...
## <20 more elements>Now we can create matrix of read enrichment at the positions of interest. We are using functions from ChIPpeakAnno package.
# load library
library(ChIPpeakAnno)
# featureAlignedSignal needs the coverage (cvg) stored in a list.
cvglist <- list(BRG1=cvg.BRG1)
sig <- featureAlignedSignal(cvglist, centers.BRG1,n.tile=200,
upstream=5000, downstream=5000)
dim(sig$BRG1)## [1] 30874 2002.7.2.2 Create heatmaps
Let’s draw the heatmaps using the EnrichedHeatmap library.
## Load the library
library(EnrichedHeatmap)
## Create a list of normalizedMatrix that is the input format
## for EnrichedHeatmap
mat1 <- list()
mat1[["BRG1"]] <- as.normalizedMatrix(as.matrix(sig[["BRG1"]]),
k_upstream = 100,
k_downstream = 100,
k_target = 0,
extend = c(5000, 5000),
signal_name = names(sig[["BRG1"]]),
target_name = "Peak centers"
)
## Create the Heatmap with default parameters
EnrichedHeatmap(mat1$BRG1, name = "BRG1")
EnrichedHeatmap combines the average profile and the density heatmap. We can observe a greater precision of the signal around the peak centers. As with ChIPpeakAnno, by default, heatmaps are sorted by read enrichment. However, it would be worth grouping together regions that have similar read enrichment pattern. This can be done using a clustering method such as k-means. This type of clustering requires the number of expected clusters to be set. Moreover, to obtain reproducible clustering results, we need to set a seed.
## define a seed value to get the same results when re-running the analysis
set.seed(123)
## Create Heatmaps with k-means clustering on BRG1 data
## We keep the generated object in order to use the clustering
## information.
heatmap.kmeans <- EnrichedHeatmap(mat1$BRG1, name = "BRG1", row_km = 8,
column_title = "BRG1", row_title_rot = 0)
## draw the heatmap
htlist <- draw(heatmap.kmeans)## Warning: did not converge in 10 iterations
Let’s enhance it!
- We first retrieve the peak clusters from
htlistobject to use them as a pre-defined partition in our enhanced heatmap.
# Use the row_order function to retrieve peaks index
# belonging to each cluster
clusters <- row_order(htlist)
# rename the clusters
names(clusters) <- paste0("cluster",names(clusters))
# transform to a vector
# Check the class and length of the resulting object
class(clusters)## [1] "list"length(clusters)## [1] 8# Check what is actually in the list elements
head(clusters[[1]])## [1] 1353 1815 22265 1564 19335 19294# Each element of the list contain the indexes of the peaks in the
# original object
# Check cluster sizes
lapply(clusters,length)## $cluster1
## [1] 1214
##
## $cluster2
## [1] 1462
##
## $cluster3
## [1] 1691
##
## $cluster4
## [1] 3055
##
## $cluster5
## [1] 2035
##
## $cluster6
## [1] 5363
##
## $cluster7
## [1] 15114
##
## $cluster8
## [1] 940# create the partition by transforming the list in vector
partition <- unlist(clusters)
head(partition)## cluster11 cluster12 cluster13 cluster14 cluster15 cluster16
## 1353 1815 22265 1564 19335 19294# names of the elements were extended with a number
# We thus trim them to retrieve the cluster name
names(partition) <- substring(first = 1,last = 8,text=names(partition))
# The numbers are the indexes of the rows, we need to
# sort the indexes to get the right order of cluster
# labels
partition <- names(partition)[order(partition)]
# keep this information with you BRG1 peaks
peaks$BRG1$cluster <- partition
head(peaks$BRG1)## GRanges object with 6 ranges and 9 metadata columns:
## seqnames ranges strand | V4 V5
## <Rle> <IRanges> <Rle> | <character> <integer>
## [1] chr1 983476-984526 * | BRG1siCTRL_CHIP-seq_.. 108
## [2] chr1 1000729-1001179 * | BRG1siCTRL_CHIP-seq_.. 90
## [3] chr1 1020914-1021207 * | BRG1siCTRL_CHIP-seq_.. 114
## [4] chr1 1069156-1070192 * | BRG1siCTRL_CHIP-seq_.. 96
## [5] chr1 1079517-1080087 * | BRG1siCTRL_CHIP-seq_.. 179
## [6] chr1 1304853-1306013 * | BRG1siCTRL_CHIP-seq_.. 108
## V6 V7 V8 V9 V10 center.peak
## <character> <numeric> <numeric> <numeric> <integer> <numeric>
## [1] . 5.97792 13.5966 10.83201 384 984001
## [2] . 6.65677 11.6289 9.01201 199 1000954
## [3] . 8.11634 14.2980 11.46427 158 1021060
## [4] . 6.90466 12.3485 9.65604 410 1069674
## [5] . 9.95583 21.2517 17.96153 213 1079802
## [6] . 5.97792 13.5966 10.83201 176 1305433
## cluster
## <character>
## [1] cluster7
## [2] cluster7
## [3] cluster7
## [4] cluster7
## [5] cluster7
## [6] cluster6
## -------
## seqinfo: 24 sequences from an unspecified genome; no seqlengths- Now, customize the enrichment colors and differentiate clusters with colors.
library(circlize)
# Define new colors for each heatmap
col_brg1 <- colorRamp2(c(0,10,15), c("white", "blue","black"))
# create a legend for the cluster labels
lgd <- Legend(at = c("cluster1", "cluster2", "cluster3", "cluster4","cluster5","cluster6","cluster7","cluster8"),
title = "Clusters",
type = "lines", legend_gp = gpar(col = 2:9))
# Add a first column containing the cluster assignment
ht_list <- Heatmap(partition, col = structure(2:9, names = paste0("cluster", 1:8)), name = "partition", show_row_names = FALSE, width = unit(3, "mm")) +
EnrichedHeatmap(mat1$BRG1, name = "BRG1", col=col_brg1, # specify the heat colors
# color per cluster
top_annotation = HeatmapAnnotation(lines = anno_enriched(gp = gpar(col = 2:9))),
column_title = "BRG1")
draw(ht_list, split = partition, annotation_legend_list = list(lgd),
ht_gap = unit(c(2, 8), "mm"))
2.7.2.3 Are BRG1 binding profiles associated with particular genomic regions ?
We will use here the clusters obtained with EnrichedHeatmap and look at the genomic distribution of peaks using TxDb.Hsapiens.UCSC.hg38.knownGene and plotAnnoBar.
# The following code uses each element of the list (peaks indexes)
# to select corresponding peaks in the original data and return them
BRG1.clusters <- lapply(clusters,function(x,peaks){
# x represents one element of the list, one set of
# indexes
cl <- peaks[x]
return(cl)
},peaks=peaks$BRG1)
# check names of the resulting list
names(BRG1.clusters)## [1] "cluster1" "cluster2" "cluster3" "cluster4" "cluster5" "cluster6" "cluster7"
## [8] "cluster8"# Transform the list as GRangesList to be able to use it with the annotatePeak
# function
BRG1.clusters <- as(BRG1.clusters,"GRangesList")
# apply annotatePeak function to each element of the BRG1.clusters GRangesList
peakAnnoList <- lapply(BRG1.clusters, annotatePeak, TxDb=txdb,
tssRegion=c(-1000, 100), annoDb="org.Hs.eg.db",
verbose=FALSE)
plotAnnoBar(peakAnnoList)
3 Data integration: let’s compare BRG1 to MITF and SOX10
3.1 Dataset description
3.1.1 Load MITF and SOX10 datasets
## Let's load MITF and SOX10 peak sets
peaks[["MITF"]] <- readPeakFile(paste0(params$chipseq, params$mitf), as="GRanges")
peaks[["SOX10"]] <- readPeakFile(paste0(params$chipseq, params$sox10), as="GRanges")3.1.2 How many peaks were called?
Compute the number of datasets.
length(peaks)## [1] 3# Load ggplot2 library
library(ggplot2)
peak.lengths <- data.frame(IP=names(peaks),
NbPeaks=sapply(peaks,length))
# check the object
peak.lengths## IP NbPeaks
## BRG1 BRG1 30874
## MITF MITF 9702
## SOX10 SOX10 4538ggplot(peak.lengths, aes(x=IP, y=NbPeaks, fill=IP)) +
geom_bar(stat="identity")+
scale_fill_brewer(palette="Set1")+
theme_classic()
3.1.3 How large are these peaks?
Compute statistics on all peak sizes
peak.width = lapply(peaks,width)
str(peak.width)## List of 3
## $ BRG1 : int [1:30874] 1051 451 294 1037 571 1161 360 285 897 905 ...
## $ MITF : int [1:9702] 260 190 262 175 195 257 233 188 1251 317 ...
## $ SOX10: int [1:4538] 252 272 331 517 253 362 246 306 254 400 ...peak.width.table <- melt(peak.width)
head(peak.width.table)## value L1
## 1 1051 BRG1
## 2 451 BRG1
## 3 294 BRG1
## 4 1037 BRG1
## 5 571 BRG1
## 6 1161 BRG1In this table, there is as many rows as the total number of peaks in all peak sets. It contains all possible pairs IP <-> number of peaks possible.
# Number of peaks per chipped protein
sapply(peaks, length)## BRG1 MITF SOX10
## 30874 9702 4538# total number of peaks
sum(sapply(peaks, length))## [1] 45114# size of the table we've just generated
dim(peak.width.table)## [1] 45114 2ggplot(peak.width.table, aes(x=L1, y=value, fill=L1)) +
geom_boxplot()+
theme_classic()+
scale_fill_brewer(palette="Set1")+
labs(x = "", y = "Peak sizes", fill = "")+
scale_y_log10()
3.2 Peak filtering
To make sure we keep only high quality data. We are going to select one those peaks having a qValue >= 4. The qValue corresponds to the 9th column of narrowPeak files. So, we are going to set a threshold on this.
## Select high quality peaks in each element of the list peaks
peaks <- lapply(peaks, function(x){
x[x$V9 >= 4,]
})
## Compute the number of remaining peaks
sapply(peaks, length)## BRG1 MITF SOX10
## 30874 9279 42143.2.1 Peak annotation
Load and add MITF and SOX10 annotations to the peakAnno list already containing BRG1 annotations.
peakAnno[["MITF"]] = annotatePeak(peaks$MITF, tssRegion=c(-1000, 100), TxDb=txdb, annoDb="org.Hs.eg.db")## >> preparing features information... 2021-05-21 19:54:40
## >> identifying nearest features... 2021-05-21 19:54:40
## >> calculating distance from peak to TSS... 2021-05-21 19:54:41
## >> assigning genomic annotation... 2021-05-21 19:54:41
## >> adding gene annotation... 2021-05-21 19:54:45
## >> assigning chromosome lengths 2021-05-21 19:54:45
## >> done... 2021-05-21 19:54:45peakAnno[["SOX10"]] = annotatePeak(peaks$SOX10, tssRegion=c(-1000, 100), TxDb=txdb, annoDb="org.Hs.eg.db")## >> preparing features information... 2021-05-21 19:54:45
## >> identifying nearest features... 2021-05-21 19:54:45
## >> calculating distance from peak to TSS... 2021-05-21 19:54:45
## >> assigning genomic annotation... 2021-05-21 19:54:45
## >> adding gene annotation... 2021-05-21 19:54:49
## >> assigning chromosome lengths 2021-05-21 19:54:49
## >> done... 2021-05-21 19:54:49plotAnnoBar(peakAnno)
3.3 Are BRG1, MITF and SOX10 co-localizing ?
3.3.1 Compare BRG1, MITF and SOX10 peak positions (Venn Diagram)
We can evaluate the peak overlapping with a Venn diagram using the ChIPpeakAnno package.
library(ChIPpeakAnno)
# We first compute the overlap between peak sets, keeping the information
# of all peaks overlapping in each set (see ?findOverlapsOfPeaks for help)
ovl <- findOverlapsOfPeaks(peaks, connectedPeaks="keepAll")ChIPpeakAnno imposes, while plotting the Venn diagram, to compute the significance of the pairwise associations using a hypergeometric test. To this end, we need to estimate the number of all potential binding events which is used by the makeVennDiagram function through the totalTest number. It is used for the hypergeometric sampling that is used to determine if the overlap between two datasets is more than would be expected by chance. This is not a trivial question, the answer is driven by what you know about the binding properties of your factors (eg. sequence specific, mainly intergenic etc). You can find an interesting discussion here. In our case we can refer to the genomic distribution of the peaks that we have plotted previously. We can assume that our TFs have a gene body binding preference. Genes cover roughly 10% of the genome.
# Estimate the average size of the peaks ...
averagePeakWidth <- mean(width(unlist(GRangesList(ovl$peaklist))))
# ... to count how many potential sites could have been bound in coding regions.
tot <- ceiling(3.3e+9 * 0.1 / averagePeakWidth)TIPS: We can define the colors attributed to each set using the function colours(). In any case you want to set a color you can use this function. Please, have a look to this page or more generally this page if you are interested in finely customizing colors of your plots.
makeVennDiagram(ovl, totalTest=tot, connectedPeaks="keepAll",
fill=brewer.pal(3,"Set1"), # circle fill color
col=brewer.pal(3,"Set1"), #circle border color
cat.col=brewer.pal(3,"Set1"))
## $p.value
## BRG1 MITF SOX10 pval
## [1,] 0 1 1 7.824528e-215
## [2,] 1 0 1 0.000000e+00
## [3,] 1 1 0 0.000000e+00
##
## $vennCounts
## BRG1 MITF SOX10 Counts count.BRG1 count.MITF count.SOX10
## [1,] 0 0 0 411865 0 0 0
## [2,] 0 0 1 1127 0 0 1127
## [3,] 0 1 0 6930 0 6930 0
## [4,] 0 1 1 40 0 40 40
## [5,] 1 0 0 24781 24781 0 0
## [6,] 1 0 1 2547 3491 0 2572
## [7,] 1 1 0 1784 1907 1822 0
## [8,] 1 1 1 454 695 487 475
## attr(,"class")
## [1] "VennCounts"According to the hypergeometric test p-values all pairwise comparisons are highly significant (=0). This has to be taken very carefully as it depends largely on the background estimation that we may have over-estimated. If you are not sure about your estimation, you will prefer to use a non-parametric approach based on your peak genomic distribution to estimate randomness. To this end the TxDb.Hsapiens.UCSC.hg38.knownGene is used.
3.3.1.1 Is the overlap significant?
ChIPpeakAnno provides the preparePool and peakPermTest functions to compute this test. These tests are made by pair, we will look at SOX10/MITF overlap significance as an example.
# Prepare a pool of random peaks following the characteristics of our peak sets
pool <- preparePool(txdb,peaks$SOX10,bindingType="TSS",featureType="transcript",seqn=paste0("chr",c(1:22,"X","Y")))
# Create the permPool object needed for peakPermTest
pool <- new("permPool",grs=pool$grs[1],N=length(peaks$SOX10))
SOX10.MITF <- peakPermTest(peaks$SOX10, peaks$MITF, pool=pool, seed=1, force.parallel=FALSE)
SOX10.MITF## $cntOverlaps
## P-value: 0.0099009900990099
## Z-score: 139.5137
## Number of iterations: 100
## Alternative: greater
## Evaluation of the original region set: 437
## Evaluation function: cntOverlaps
## Randomization function: randPeaks
##
## attr(,"class")
## [1] "permTestResultsList"plot(SOX10.MITF)
Venn diagrams are widely used to represent overlaps, intersections. However, in ChIP-seq analysis, the definition of the peaks is dependent on the peak caller, the FDR thresholds (what about peaks just bellow the threshold). The overlap is also difficult to assess, indeed do we call an overlap a 1 nucleotide an intersection ? It is one of the numerous parameters that can be tuned …
Are different combinations of TFs bind specific genomic regions ?
coocs <- as(ovl$peaklist,"GRangesList")
# apply annotatePeak to each set of peaks in the list coocs
peakAnnoList <- lapply(coocs, annotatePeak, TxDb=txdb,
tssRegion=c(-1000, 100), annoDb="org.Hs.eg.db",
verbose=FALSE)## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns# Plot peak distribution relatively to gene features
plotAnnoBar(peakAnnoList)
# Plot peak distribution relatively to their distance to the TSS
plotDistToTSS(peakAnnoList)
3.3.2 Heatmaps / Profiles
One way to circumvent hard thresholds and relatively arbitrary choices, we can choose to use heatmap and average profile representations. To this end, we need to define the reference point of view. We use the BRG1 peak centers.
Load bigwig file for MITF and SOX10 at BRG1 peaks.
# Load and add bigwig profiles for MITF and SOX10 in cvglist already containing
# BRG1
names(cvglist)## [1] "BRG1"cvglist$MITF <- import(file.path(params$chipseq,params$mitf.bw),format="BigWig",
which=extended.5K.BRG1,
as="RleList")
cvglist$SOX10 <- import(file.path(params$chipseq,params$sox10.bw),format="BigWig",
which=extended.5K.BRG1,
as="RleList")
names(cvglist)## [1] "BRG1" "MITF" "SOX10"Prepare the matrices binned in 50bp windows.
sig <- featureAlignedSignal(cvglist, centers.BRG1,
upstream=5000, downstream=5000,n.tile=200)
lapply(sig, dim)## $BRG1
## [1] 30874 200
##
## $MITF
## [1] 30874 200
##
## $SOX10
## [1] 30874 200Create normalized matrices.
## Create a list of normalizedMatrix that is the input format
## for EnrichedHeatmap
BRG1.mat <- lapply(sig, function(x){
# x represent each element of the sig list
mat <- as.normalizedMatrix(as.matrix(x),
k_upstream = 100,
k_downstream = 100,
k_target = 0,
extend = c(5000, 5000),
#signal_name = names(sig[["MITF"]]),
target_name = "Peak center"
)
return(mat)
})We use the partition we computed in section 2.7.2.2
# Define new colors for each heatmap
col_sox10 <-colorRamp2(c(0,3,4), c("white", "blue","black"))
col_mitf <- colorRamp2(c(0,4,5), c("white", "blue","black"))
col_brg1 <- colorRamp2(c(0,5,6), c("white", "blue","black"))
# create a legend for the cluster labels
lgd <- Legend(at = c("cluster1", "cluster2", "cluster3", "cluster4","cluster5","cluster6","cluster7","cluster8"),
title = "Clusters",
type = "lines", legend_gp = gpar(col = 2:9))
# Add a first column containing the cluster assignment
ht_list <- Heatmap(partition, col = structure(2:9, names = paste0("cluster", 1:8)), name = "partition", show_row_names = FALSE, width = unit(3, "mm")) +
EnrichedHeatmap(BRG1.mat$BRG1, name = "BRG1", col=col_brg1,
top_annotation = HeatmapAnnotation(lines = anno_enriched(gp = gpar(col = 2:9))),
column_title = "BRG1") +
EnrichedHeatmap(BRG1.mat$MITF, name = "MITF",
top_annotation = HeatmapAnnotation(lines = anno_enriched(gp = gpar(col = 2:9))),
column_title = "MITF", row_title_rot = 0, col=col_mitf) +
EnrichedHeatmap(BRG1.mat$SOX10, name = "SOX10",
top_annotation = HeatmapAnnotation(lines = anno_enriched(gp = gpar(col = 2:9))),
column_title = "SOX10", row_title_rot = 0, col=col_sox10)
draw(ht_list, split = partition, annotation_legend_list = list(lgd),
ht_gap = unit(c(2, 8, 8,8), "mm"))
Lets have a look at the TF binding profiles in cluster 5.
# For each matrix, filter peaks from cluster 5
sig.tmp <- lapply(sig,function(x,part){
mat <- x[part=="cluster5",]
return(mat)
},part=partition)
# plot the average profiles
featureAlignedDistribution(sig.tmp, centers.BRG1,
upstream=5000, downstream=5000,
type="l",n.tile=200)
3.4 Is gene expression influenced by TF binding to their promoters ?
- Select genes having BRG1 binding close to their TSS. To do that we will plot a heatmap centered on TSS.
First retrieve the TSS positions from the TxDB object
# Generate TSS and promoter GRanges, the function promoters allows to retrieve the
# gene ID (entrez)
TSS <- GenomicFeatures::promoters(txdb, upstream=0, downstream=0, columns="gene_id")Clean the TSS list removing TSS with no gene ID assigned (column gene_id in the metadata) and on scaffolds (seqnames).
The column gene_id is of type characterList, this type is difficult to manipulate. We then decide to transform it to a character vector. We can’t directly unlist the object because, to save memory, not annotated TSS are assigned empty elements which disappear when unlisting.
# For each element of the characterList we check it's length, if it is 0
# we fill the element with NA. We the unlist.
gene_ids <- unlist(lapply(TSS$gene_id,function(x){
if (length(x)==0){
NA
}
else {
x
}
} ))
# Keep TSS on selected chromosomes annotated with a gene ID and remove
# duplicate IDs due to transcript isoforms. Here, the selection of the
# isoform is random.
TSS <- TSS[as.vector(seqnames(TSS))%in%paste0("chr",c(1:22,"X","Y")) & !is.na(gene_ids) & !duplicated(gene_ids)]
# extend TSS both direction
TSS.extended <- resize(TSS, width = width(TSS)+4000, fix = "center")Load bigwig profiles for TSS extended regions, create matrices and transform in the enrichedHeatmap recognized format.
# load the bw file for all TFs within the promoter
cvglist <- list()
cvglist$BRG1 <- import(file.path(params$chipseq,params$brg1.bw),format="BigWig",
which=TSS.extended,
as="RleList")
cvglist$MITF <- import(file.path(params$chipseq,params$mitf.bw),format="BigWig",
which=TSS.extended,
as="RleList")
cvglist$SOX10 <- import(file.path(params$chipseq,params$sox10.bw),format="BigWig",
which=TSS.extended,
as="RleList")
# Produce the signal matrices
sig <- featureAlignedSignal(cvglist, TSS,n.tile=400,
upstream=2000, downstream=2000)
# Transform the signal matrices as normalizedMatrix
TSS.mat <- lapply(sig, function(x){
# x represent each element of the sig list
mat <- as.normalizedMatrix(as.matrix(x),
k_upstream = 200,
k_downstream = 200,
k_target = 0,
extend = c(2000, 2000),
target_name = "TSS"
)
return(mat)
})We will compute K-means clustering with kmeans function to define the partitions of the TSS based on BRG1 signal.
# Compute a partition using the kmeans function, asking for 5 clusters
set.seed(20210526)
partition.TSS = paste0("cluster", kmeans(TSS.mat$BRG1, centers = 5)$cluster)
# Specify colors for the position enrichment for each matrix
col_sox10 = colorRamp2(c(0,3,4), c("white", "blue","black"))
col_mitf = colorRamp2(c(0,4,5), c("white", "blue","black"))
col_brg1 = colorRamp2(c(0,5,6), c("white", "blue","black"))
lgd = Legend(at = c("cluster1", "cluster2", "cluster3", "cluster4","cluster5"),
title = "Clusters",
type = "lines", legend_gp = gpar(col = 2:6))
ht_list = Heatmap(partition.TSS, col = structure(2:7, names = paste0("cluster", 1:5)), name = "partition", show_row_names = FALSE, width = unit(3, "mm")) +
EnrichedHeatmap(TSS.mat$BRG1, name = "BRG1", col=col_brg1,
top_annotation = HeatmapAnnotation(lines = anno_enriched(gp = gpar(col = 2:6))),
column_title = "BRG1") +
EnrichedHeatmap(TSS.mat$MITF, name = "MITF",
top_annotation = HeatmapAnnotation(lines = anno_enriched(gp = gpar(col = 2:6))),
column_title = "MITF", row_title_rot = 0, col=col_mitf) +
EnrichedHeatmap(TSS.mat$SOX10, name = "SOX10",
top_annotation = HeatmapAnnotation(lines = anno_enriched(gp = gpar(col = 2:6))),
column_title = "SOX10", row_title_rot = 0, col=col_sox10)
draw(ht_list, split = partition.TSS, annotation_legend_list = list(lgd),
ht_gap = unit(c(2, 8, 8,8), "mm"))
Lets have a look at the TF binding profiles in cluster 4.
# For each matrix, filter peaks from cluster 5
sig.tmp <- lapply(sig,function(x,part){
mat <- x[part=="cluster4",]
return(mat)
},part=partition.TSS)
# plot the average profiles
featureAlignedDistribution(sig.tmp, TSS,
upstream=2000, downstream=2000,
type="l",n.tile=400)## Warning in featureAlignedDistribution(sig.tmp, TSS, upstream = 2000, downstream
## = 2000, : feature.gr is set to the center of feature.gr
- Lets’ integrate expression data from RNA-seq experiment.
The goal is to visualize the gene expression according to the cluster the TSS belong.
# load RNA-seq normalized expression data
load(params$rnaseq) # RNAseq
head(RNAseq)## Ensembl.Gene.ID WT.norm shBGR1.norm logCPM logFC FDR
## 1 ENSG00000000003 722.38046 119.10483 4.387012 2.5748948 2.233826e-21
## 2 ENSG00000000419 1481.30280 2514.73393 6.689308 -0.6994493 1.391964e-13
## 3 ENSG00000000457 259.21973 314.27718 5.364018 -0.3167911 3.091766e-05
## 4 ENSG00000000460 333.22470 380.48495 5.562143 -0.1954216 2.875151e-03
## 5 ENSG00000000971 13.42921 13.59959 0.341927 0.6570620 2.536438e-02
## 6 ENSG00000001036 1157.52492 1521.75139 5.567526 -0.3587066 3.117933e-06
## signif
## 1 stable
## 2 stable
## 3 stable
## 4 stable
## 5 stable
## 6 stableRNA-seq columns: - Ensembl.Gene.ID: Gene idenfier from ensEMBL (ENSG*****) - WT.norm: normalized expression for wild type cell lines - shBGR1.norm: normalized expression for shBRG1 cell lines - logCPM: average log transformed overall expression - logFC: log fol change between the two conditions
- FDR: multi-testing corrected p-value from edgeR statistical test - signif (up, down, stable): significance was set as: FDR < 0.01, logCPM > 5, abs(logFC) > 1
Reads from RNA-seq were assigned to genes modeled from ensEMBL and thus are identified with ensEMBL ID. However, the TSS are annotated with ENTREZ IDs, we thus need to map together the identifiers. We use an OrgDB object for the human genome the available information are displayed using columns().
columns(org.Hs.eg.db)## [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
## [6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
## [11] "GO" "GOALL" "IPI" "MAP" "OMIM"
## [16] "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM" "PMID"
## [21] "PROSITE" "REFSEQ" "SYMBOL" "UCSCKG" "UNIGENE"
## [26] "UNIPROT"We then use the function mapIds from AnnotationDbi package with the keys ENSEMBL and ENTREZID. You can see that other identifiers are available from this object (e.g: SYMBOL, REFSEQ, UNIPROT etc).
conv <- mapIds(x=org.Hs.eg.db,
keys=RNAseq$Ensembl.Gene.ID, # what do we want to be mapped
column="ENTREZID", # which type of ID we want
keytype="ENSEMBL") # what type of ID we give## 'select()' returned 1:many mapping between keys and columnshead(conv)## ENSG00000000003 ENSG00000000419 ENSG00000000457 ENSG00000000460 ENSG00000000971
## "7105" "8813" "57147" "55732" "3075"
## ENSG00000001036
## "2519"# names of the vector are ENSEMBL IDs and elements are ENTREZ IDsWe use the match function that allows to return match between the IDs in RNAseq and the ensEMBL IDs from the conversion vector
m <- match(RNAseq$Ensembl.Gene.ID,names(conv))
head(m)## [1] 1 2 3 4 5 6# m contains NA for element in the first vector not present in the 2nd
# or the index of the element of the 2nd vector corresponding the those
# in the first
# Attribute the matched ENTREZ ID
RNAseq$ENTREZID <- conv[m]
head(RNAseq)## Ensembl.Gene.ID WT.norm shBGR1.norm logCPM logFC FDR
## 1 ENSG00000000003 722.38046 119.10483 4.387012 2.5748948 2.233826e-21
## 2 ENSG00000000419 1481.30280 2514.73393 6.689308 -0.6994493 1.391964e-13
## 3 ENSG00000000457 259.21973 314.27718 5.364018 -0.3167911 3.091766e-05
## 4 ENSG00000000460 333.22470 380.48495 5.562143 -0.1954216 2.875151e-03
## 5 ENSG00000000971 13.42921 13.59959 0.341927 0.6570620 2.536438e-02
## 6 ENSG00000001036 1157.52492 1521.75139 5.567526 -0.3587066 3.117933e-06
## signif ENTREZID
## 1 stable 7105
## 2 stable 8813
## 3 stable 57147
## 4 stable 55732
## 5 stable 3075
## 6 stable 2519We can now associate RNA-seq and TSS information using the ENTREZ IDs and the match function.
# We add the partition (clusters) to the TSS object in order to add the
# clusters to the RNA-seq object
TSS$cluster <- partition.TSS
# The common IDs this time are the ENTREZ IDs
m <- match(RNAseq$ENTREZID,unlist(TSS$gene_id))
RNAseq$cluster <- TSS$cluster[m]
head(RNAseq)## Ensembl.Gene.ID WT.norm shBGR1.norm logCPM logFC FDR
## 1 ENSG00000000003 722.38046 119.10483 4.387012 2.5748948 2.233826e-21
## 2 ENSG00000000419 1481.30280 2514.73393 6.689308 -0.6994493 1.391964e-13
## 3 ENSG00000000457 259.21973 314.27718 5.364018 -0.3167911 3.091766e-05
## 4 ENSG00000000460 333.22470 380.48495 5.562143 -0.1954216 2.875151e-03
## 5 ENSG00000000971 13.42921 13.59959 0.341927 0.6570620 2.536438e-02
## 6 ENSG00000001036 1157.52492 1521.75139 5.567526 -0.3587066 3.117933e-06
## signif ENTREZID cluster
## 1 stable 7105 cluster3
## 2 stable 8813 cluster4
## 3 stable 57147 cluster4
## 4 stable 55732 cluster3
## 5 stable 3075 cluster3
## 6 stable 2519 cluster4We first want to compare the expression between this to condition in the different clusters.
# melt the RNAseq table to get a data.frame with a column indicating
# the experimental condition, a column with the cluster assignment
# and the expression values
expr.melt <- melt(RNAseq[,c("WT.norm","shBGR1.norm","cluster")])## Using cluster as id variableshead(expr.melt)## cluster variable value
## 1 cluster3 WT.norm 722.38046
## 2 cluster4 WT.norm 1481.30280
## 3 cluster4 WT.norm 259.21973
## 4 cluster3 WT.norm 333.22470
## 5 cluster3 WT.norm 13.42921
## 6 cluster4 WT.norm 1157.52492ggplot(expr.melt,aes(x=cluster,y=log2(value+1))) +
geom_boxplot(aes(color=variable)) +
theme_bw() + theme(legend.position = "top") +
scale_color_manual(values=c(WT.norm="black",shBGR1.norm=colours()[613]))
We can see that there are not apparent difference between the conditions even if cluster 3, which has the less BRG1 signal, is much less expressed than the other. We can check if the distribution of significantly differentially expressed genes differs among the clusters.
# Distribution of differentially expressed genes as barplots
ggplot(RNAseq,aes(x=cluster)) +
geom_bar(aes(fill=signif),position = "fill",color="black") +
theme_bw() +
scale_fill_manual(values=c(up="red",down="blue",stable="white"))
3.5 Are the differentially expressed genes with different BRG1 binding patterns have different biological function ?
To answer this question, we use the R package clusterProfiler to compute hypergeometric enrichment of biological function in each cluster using KEGG annotation database.
clusterProfiler provide the function compareCluster that allows to analyze and compare enrichment in different group of genes. This function recognize ENTREZ identifiers.
Représentation simplifiée graphique du test hypergéometrique
# load the library
library(clusterProfiler)
# Select genes whose gene expression is significantly changing
# between conditions
geneList <- RNAseq[RNAseq$signif!="stable",]
# The function split help at splitting a vector to a list following categories
# in an other vector.
geneList <- split(geneList$ENTREZID,geneList$cluster)
names(geneList)## [1] "cluster1" "cluster2" "cluster3" "cluster4" "cluster5"head(geneList$cluster1)## [1] "27072" "54443" "57514" "57609" "998" "58476"compKEGG <- compareCluster(geneCluster = geneList,
fun = "enrichKEGG",
pvalueCutoff = 0.05,
pAdjustMethod = "BH")
dotplot(compKEGG, showCategory = 15, title = "KEGG Pathway Enrichment Analysis")
Let’s focus on cluster 4 and visualize which genes of the the KEGG pathway “Cell cycle” (hsa04110) are both differentially expressed and bound by BRG1 and MITF. We use the pathview package which maps differential gene expression values to KEGG maps and create a PNG file in your current directory that you get with getwd().
# load the library
library(pathview)
# Retrieve gene id from cluster 4
geneList.cl4 <- geneList$cluster4
# use match function to select correspondant logFC from RNA-seq data
m <- match(geneList.cl4,RNAseq$ENTREZID)
logFC.cl4 <- RNAseq$logFC[m]
names(logFC.cl4) <- geneList.cl4
# call the function
pathview(logFC.cl4,pathway.id = "hsa04110",species = "hsa")
3.5.1 Are there published data I could use to enrich my analysis ?
- Data Mining with ChIP seq data deposited in GEO with
ChIPSeeker
There are many ChIP seq data sets that have been published and deposited in GEO database. We can compare our own dataset to those deposited in GEO to search for significant overlap data. Significant overlap of ChIP seq data by different binding proteins may be used to infer cooperative regulation and thus can be used to generate hypotheses.
We collect about 17,000 bed files deposited in GEO, user can use getGEOspecies to get a summary based on spieces.
getGEOspecies()## species Freq
## 1 Aedes aegypti 11
## 2 Anabaena 6
## 3 Anolis carolinensis 5
## 4 Anopheles gambiae 2
## 5 Apis mellifera 5
## 6 Apis mellifera scutellata 1
## 7 Arabidopsis lyrata 4
## 8 Arabidopsis thaliana 288
## 9 Atelerix albiventris 2
## 10 Bos taurus 37
## 11 Branchiostoma lanceolatum 62
## 12 Brassica rapa 12
## 13 Caenorhabditis elegans 189
## 14 Candida albicans 25
## 15 Candida dubliniensis 20
## 16 Canis lupus familiaris 7
## 17 Chlamydomonas reinhardtii 51
## 18 Chlorocebus aethiops 2
## 19 Cleome hassleriana 1
## 20 Columba livia 6
## 21 Crassostrea gigas 1
## 22 Cryptococcus neoformans 51
## 23 Cyprinus carpio 40
## 24 Danio rerio 308
## 25 Drosophila busckii 1
## 26 Drosophila melanogaster 1069
## 27 Drosophila melanogaster;\tSindbis virus 3
## 28 Drosophila miranda 2
## 29 Drosophila pseudoobscura 7
## 30 Drosophila simulans 12
## 31 Drosophila virilis 26
## 32 Drosophila willistoni 1
## 33 Drosophila yakuba 8
## 34 Equus caballus 1
## 35 Escherichia coli 15
## 36 Escherichia coli BW25113 4
## 37 Escherichia coli K-12 2
## 38 Escherichia coli str. K-12 substr. MG1655 16
## 39 Gallus gallus 58
## 40 Geobacter sulfurreducens PCA 3
## 41 Gorilla gorilla 2
## 42 Histophilus somni 1
## 43 Homo sapiens 29978
## 44 Homo sapiens;\tHuman herpesvirus 8 6
## 45 Human herpesvirus 6B 2
## 46 Human herpesvirus 8 6
## 47 Larimichthys crocea 7
## 48 Legionella pneumophila 5
## 49 Leishmania amazonensis 4
## 50 Leishmania major 2
## 51 Leishmania major strain Friedlin 4
## 52 Leishmania tarentolae 15
## 53 Macaca mulatta 120
## 54 Monodelphis domestica 8
## 55 Moraxella catarrhalis O35E 6
## 56 Mus 2
## 57 Mus musculus 16748
## 58 Mus musculus x Mus spretus 1
## 59 Mycobacterium tuberculosis 2
## 60 Myotis brandtii 15
## 61 Naumovozyma castellii 1
## 62 Nematostella vectensis 23
## 63 Oreochromis niloticus 1
## 64 Ornithorhynchus anatinus 5
## 65 Oryza sativa 30
## 66 Oryzias latipes 2
## 67 Pan troglodytes 93
## 68 Papio anubis 1
## 69 Plasmodium falciparum 129
## 70 Plasmodium falciparum 3D7 29
## 71 Pseudomonas putida KT2440 2
## 72 Pseudozyma aphidis 11
## 73 Pyrococcus furiosus 4
## 74 Rattus norvegicus 108
## 75 Rhodopseudomonas palustris 6
## 76 Rhodopseudomonas palustris CGA009 3
## 77 Saccharomyces cerevisiae 813
## 78 Saccharomyces cerevisiae x Saccharomyces paradoxus 16
## 79 Saccharomyces cerevisiae;\tMus musculus 12
## 80 Saccharomyces kudriavzevii 1
## 81 Saccharomyces paradoxus 8
## 82 Saccharomyces uvarum 1
## 83 Schizosaccharomyces japonicus 2
## 84 Schizosaccharomyces pombe 179
## 85 Schmidtea mediterranea 7
## 86 Solanum lycopersicum 2
## 87 Sorghum bicolor 2
## 88 Spodoptera frugiperda 16
## 89 Streptomyces coelicolor A3(2) 6
## 90 Sus scrofa 41
## 91 Taeniopygia guttata 1
## 92 Tupaia chinensis 7
## 93 Vibrio cholerae 8
## 94 Xenopus (Silurana) tropicalis 62
## 95 Xenopus laevis 10
## 96 Xenopus tropicalis 74
## 97 Zea mays 65User can access the detail information by getGEOInfo, for each genome version.
hg38 <- getGEOInfo(genome="hg38", simplify=TRUE)
head(hg38)## series_id gsm organism title
## 16488 GSE58207 GSM1403308 Homo sapiens Lactimidomycin treated HCT116
## 16489 GSE58207 GSM1403308 Homo sapiens Lactimidomycin treated HCT116
## 16490 GSE58207 GSM1403307 Homo sapiens Cycloheximide treated HCT116
## 16491 GSE58207 GSM1403307 Homo sapiens Cycloheximide treated HCT116
## 20345 GSE67978 GSM1660032 Homo sapiens H3K27ac_Human_Brain_WhiteMatter_HS2
## 20351 GSE67978 GSM1660029 Homo sapiens H3K27ac_Human_Brain_OccipitalPole_HS2
## supplementary_file
## 16488 ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1403nnn/GSM1403308/suppl/GSM1403308_HCT116_LTM_sense.bedGraph.gz
## 16489 ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1403nnn/GSM1403308/suppl/GSM1403308_HCT116_LTM_anti-sense.bedGraph.gz
## 16490 ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1403nnn/GSM1403307/suppl/GSM1403307_HCT116_CHX_anti-sense.bedGraph.gz
## 16491 ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1403nnn/GSM1403307/suppl/GSM1403307_HCT116_CHX_sense.bedGraph.gz
## 20345 ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1660nnn/GSM1660032/suppl/GSM1660032_H3K27ac_Human_Brain_WhiteMatter_HS2_hg38_peaks.narrowPeak.gz
## 20351 ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM1660nnn/GSM1660029/suppl/GSM1660029_H3K27ac_Human_Brain_OccipitalPole_HS2_hg38_peaks.narrowPeak.gz
## genomeVersion pubmed_id
## 16488 hg38 <NA>
## 16489 hg38 <NA>
## 16490 hg38 <NA>
## 16491 hg38 <NA>
## 20345 hg38 <NA>
## 20351 hg38 <NA>ChIPseeker provide function downloadGEObedFiles to download all the bed files of a particular genome.
downloadGEObedFiles(genome="hg38", destDir="hg38")Or a vector of GSM accession number by downloadGSMbedFiles.
gsm <- hg38$gsm[sample(nrow(hg38), 10)]
downloadGSMbedFiles(gsm, destDir="hg38")After download the bed files from GEO, we can pass them to enrichPeakOverlap for testing the significant of overlap. Parameter targetPeak can be the folder, e.g. hg19, that containing bed files. enrichPeakOverlap will parse the folder and compare all the bed files. It is possible to test the overlap with bed files that are mapping to different genome or different genome versions, enrichPeakOverlap provide a parameter chainFile that can pass a chain file and liftOver the targetPeak to the genome version consistent with queryPeak. Signifcant overlap can be use to generate hypothesis of cooperative regulation.By mining the data deposited in GEO, we can identify some putative complex or interacted regulators in gene expression regulation or chromosome remodelling for further validation.
4 Session info
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
##
## attached base packages:
## [1] parallel stats4 grid stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] pathview_1.30.1
## [2] clusterProfiler_3.18.0
## [3] UpSetR_1.4.0
## [4] reshape2_1.4.4
## [5] RColorBrewer_1.1-2
## [6] knitr_1.30
## [7] org.Hs.eg.db_3.12.0
## [8] TxDb.Hsapiens.UCSC.hg38.knownGene_3.10.0
## [9] GenomicFeatures_1.42.1
## [10] AnnotationDbi_1.52.0
## [11] Biobase_2.50.0
## [12] rtracklayer_1.50.0
## [13] circlize_0.4.12
## [14] ChIPseeker_1.26.2
## [15] ChIPpeakAnno_3.24.2
## [16] ggplot2_3.3.3
## [17] EnrichedHeatmap_1.20.0
## [18] GenomicRanges_1.42.0
## [19] GenomeInfoDb_1.26.2
## [20] IRanges_2.24.1
## [21] S4Vectors_0.28.1
## [22] BiocGenerics_0.36.0
## [23] ComplexHeatmap_2.6.2
##
## loaded via a namespace (and not attached):
## [1] shadowtext_0.0.7
## [2] fastmatch_1.1-0
## [3] BiocFileCache_1.14.0
## [4] plyr_1.8.6
## [5] igraph_1.2.6
## [6] lazyeval_0.2.2
## [7] splines_4.0.3
## [8] BiocParallel_1.24.1
## [9] digest_0.6.27
## [10] ensembldb_2.14.1
## [11] htmltools_0.5.1.1
## [12] GOSemSim_2.16.1
## [13] viridis_0.5.1
## [14] GO.db_3.12.1
## [15] magrittr_2.0.1
## [16] memoise_1.1.0
## [17] BSgenome_1.58.0
## [18] cluster_2.1.0
## [19] Biostrings_2.58.0
## [20] graphlayouts_0.7.1
## [21] matrixStats_0.57.0
## [22] askpass_1.1
## [23] enrichplot_1.10.1
## [24] prettyunits_1.1.1
## [25] colorspace_2.0-0
## [26] blob_1.2.1
## [27] rappdirs_0.3.1
## [28] ggrepel_0.9.0
## [29] xfun_0.19
## [30] dplyr_1.0.2
## [31] crayon_1.3.4
## [32] RCurl_1.98-1.2
## [33] scatterpie_0.1.5
## [34] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
## [35] graph_1.68.0
## [36] survival_3.2-7
## [37] glue_1.4.2
## [38] polyclip_1.10-0
## [39] gtable_0.3.0
## [40] zlibbioc_1.36.0
## [41] XVector_0.30.0
## [42] GetoptLong_1.0.5
## [43] DelayedArray_0.16.0
## [44] Rgraphviz_2.34.0
## [45] shape_1.4.5
## [46] scales_1.1.1
## [47] DOSE_3.16.0
## [48] futile.options_1.0.1
## [49] DBI_1.1.0
## [50] Rcpp_1.0.5
## [51] plotrix_3.8-1
## [52] viridisLite_0.3.0
## [53] progress_1.2.2
## [54] clue_0.3-59
## [55] bit_4.0.4
## [56] httr_1.4.2
## [57] fgsea_1.16.0
## [58] gplots_3.1.1
## [59] ellipsis_0.3.1
## [60] pkgconfig_2.0.3
## [61] XML_3.99-0.5
## [62] farver_2.0.3
## [63] dbplyr_2.0.0
## [64] locfit_1.5-9.4
## [65] labeling_0.4.2
## [66] tidyselect_1.1.0
## [67] rlang_0.4.10
## [68] munsell_0.5.0
## [69] tools_4.0.3
## [70] downloader_0.4
## [71] generics_0.1.0
## [72] RSQLite_2.2.1
## [73] evaluate_0.14
## [74] stringr_1.4.0
## [75] yaml_2.2.1
## [76] bit64_4.0.5
## [77] tidygraph_1.2.0
## [78] caTools_1.18.2
## [79] purrr_0.3.4
## [80] KEGGREST_1.30.1
## [81] AnnotationFilter_1.14.0
## [82] ggraph_2.0.4
## [83] RBGL_1.66.0
## [84] formatR_1.7
## [85] KEGGgraph_1.50.0
## [86] DO.db_2.9
## [87] xml2_1.3.2
## [88] biomaRt_2.46.0
## [89] compiler_4.0.3
## [90] curl_4.3
## [91] png_0.1-7
## [92] tibble_3.0.4
## [93] tweenr_1.0.1
## [94] stringi_1.5.3
## [95] highr_0.8
## [96] futile.logger_1.4.3
## [97] lattice_0.20-41
## [98] ProtGenerics_1.22.0
## [99] Matrix_1.3-0
## [100] multtest_2.46.0
## [101] vctrs_0.3.6
## [102] pillar_1.4.7
## [103] lifecycle_1.0.0
## [104] BiocManager_1.30.10
## [105] GlobalOptions_0.1.2
## [106] cowplot_1.1.1
## [107] data.table_1.13.6
## [108] bitops_1.0-6
## [109] qvalue_2.22.0
## [110] R6_2.5.0
## [111] KernSmooth_2.23-18
## [112] gridExtra_2.3
## [113] codetools_0.2-18
## [114] lambda.r_1.2.4
## [115] gtools_3.8.2
## [116] boot_1.3-25
## [117] MASS_7.3-53
## [118] assertthat_0.2.1
## [119] SummarizedExperiment_1.20.0
## [120] openssl_1.4.3
## [121] rjson_0.2.20
## [122] withr_2.3.0
## [123] regioneR_1.22.0
## [124] GenomicAlignments_1.26.0
## [125] Rsamtools_2.6.0
## [126] GenomeInfoDbData_1.2.4
## [127] hms_0.5.3
## [128] ggupset_0.3.0
## [129] VennDiagram_1.6.20
## [130] tidyr_1.1.2
## [131] rvcheck_0.1.8
## [132] rmarkdown_2.6
## [133] MatrixGenerics_1.2.0
## [134] Cairo_1.5-12.2
## [135] ggforce_0.3.2References
“(2014) Mod/Mouse/humanENCODE: Blacklisted Genomic Regions for Functional Genomics Analysis - Anshul Kundaje.” n.d. Accessed August 26, 2016. https://sites.google.com/site/anshulkundaje/projects/blacklists.
Anders, Simon, Paul Theodor Pyl, and Wolfgang Huber. 2015. “HTSeq—a Python Framework to Work with High-Throughput Sequencing Data.” Bioinformatics 31 (2): 166–69. https://doi.org/10.1093/bioinformatics/btu638.
Anders, and Huber. 2010. “Differential Expression Analysis for Sequence Count Data.” Genome Biology 11.
Dobin, Alexander, Carrie A. Davis, Felix Schlesinger, Jorg Drenkow, Chris Zaleski, Sonali Jha, Philippe Batut, Mark Chaisson, and Thomas R. Gingeras. 2013. “STAR: Ultrafast Universal RNA-Seq Aligner.” Bioinformatics 29 (1): 15–21. https://doi.org/10.1093/bioinformatics/bts635.
Heinz, Sven, Christopher Benner, Nathanael Spann, Eric Bertolino, Yin C. Lin, Peter Laslo, Jason X. Cheng, Cornelis Murre, Harinder Singh, and Christopher K. Glass. 2010. “Simple Combinations of Lineage-Determining Transcription Factors Prime Cis-Regulatory Elements Required for Macrophage and B Cell Identities.” Molecular Cell 38 (4): 576–89. https://doi.org/10.1016/j.molcel.2010.05.004.
Langmead, Ben, Cole Trapnell, Mihai Pop, and Steven L. Salzberg. 2009. “Ultrafast and Memory-Efficient Alignment of Short DNA Sequences to the Human Genome.” Genome Biology 10 (3): R25. https://doi.org/10.1186/gb-2009-10-3-r25.
Laurette, Patrick, Thomas Strub, Dana Koludrovic, Céline Keime, Stéphanie Le Gras, Hannah Seberg, Eric Van Otterloo, et al. 2015. “Transcription Factor MITF and Remodeller BRG1 Define Chromatin Organisation at Regulatory Elements in Melanoma Cells.” Edited by Michael R Green. eLife 4 (March): e06857. https://doi.org/10.7554/eLife.06857.
Martin, Marcel. 2011. “Cutadapt Removes Adapter Sequences from High-Throughput Sequencing Reads.” EMBnet.journal 17 (1): pp. 10–12. http://journal.embnet.org/index.php/embnetjournal/article/view/200.
Quinlan, Aaron R., and Ira M. Hall. 2010. “BEDTools: A Flexible Suite of Utilities for Comparing Genomic Features.” Bioinformatics 26 (6): 841–42. https://doi.org/10.1093/bioinformatics/btq033.
Robinson, Mark D, Davis J McCarthy, and Gordon K Smyth. 2010. “EdgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics 26 (1): 139–40. https://doi.org/10.1093/bioinformatics/btp616.